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In the development of multiscale biological models it is crucial to establish a connection between 
discrete microscopic or mesoscopic stochastic models and macroscopic continuous descriptions based 
on cellular density. In this paper a continuous limit of a two-dimensional Cellular Potts Model 
(CPM) with excluded volume is derived, describing cells moving in a medium and reacting to each 
other through both direct contact and long range chemotaxis. The continuous macroscopic model is 
obtained as a Fokker-Planck equation describing evolution of the cell probability density function. 
All coefficients of the general macroscopic model are derived from parameters of the CPM and a very 
good agreement is demonstrated between CPM Monte Carlo simulations and numerical solution 
of the macroscopic model. It is also shown that in the absence of contact cell-cell interactions, 
the obtained model reduces to the classical macroscopic Keller-Segel model. General multiscale 
approach is demonstrated by simulating spongy bone formation from loosely packed mesenchyme 
via the intramembranous route suggesting that self-organizing physical mechanisms can account for 
this developmental process. 

PACS numbers: 87.18.Ed, 05.40.Ca, 05.65.+b, 87.18.Hf, 87.18.Bb; 87.1 



A large literature exists studying continuous limits of 
point-wise discrete microscopic models for biological sys- 
tems. For example, the classic Keller-Segel PDE model 
of chemotaxis [1] was derived from a discrete model with 
point-wise cells undergoing random walk [2-5]. How- 
ever, many biological phenomena require taking into ac- 
count the finite size of biological cells, and much less 
work has been done on deriving macroscopic limits of mi- 
croscopic models which treat cells as extended objects. 
The mesoscopic Cellular Potts Model (CPM), first in- 
troduced by Glazier and Graner [6, 7], has been used 
as a component of multiscale, experimentally motivated 
hybrid approaches, combining discrete and macroscopic 
continuous representations, to simulate, among others, 
morphological phenomena in the cellular slime mold Dic- 
tyostelium discoideum vascular development Q and 
the proximo-distal increase in the number of skeletal el- 
ements in the developing avian limb [xoj | . 

One of the earliest attempts at combining mesoscopic 
and macroscopic levels of description of cellular dynam- 
ics was described in ll| where the diffusion coefficient 
for a collection of noninteracting randomly moving cells 
was derived from a one-dimensional CPM. Recently a 
microscopic limit of subcellular elements model [l2| was 
derived in the form of an advection-diffusion partial dif- 
ferential equation for cellular density. In previous pa- 
pers HI 14 1 we studied the continuous limit of ID and 
2D models of individual cell motion in a medium, in the 
presence of an external field but without contact cell-cell 
interactions. 

This paper describes a theoretical analysis leading to 
a continuous macroscopic limit of the two-dimensional 
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FIG. 1: Cell representation in the two dimensional CPM. In 
this picture grey and white colors are used to indicate the 
cell body and ECM respectively. Cell can grow or shrink in x 
and y direction by adding or removing one row (or column) 
of pixels. 



mesocopic CPM with contact cell-cell interactions. Our 
approach, which is based on combining mesoscopic and 
macroscopic models, can be applied to studying biolog- 
ical phenomena in which a nonconfluent population of 
cells interact directly and via soluble factors, forming an 
open network structure. Examples include vasculogene- 
sis d 
bone 



, 17 1 and formation of trabecular, or spongy, 
13 . [2f| to be described below. 
The CPM, defined on a multidimensional lattice, al- 
lows simulation of both cell-cell contact and chemotactic 
long distance interactions, along with extended cell rep- 
resentations. In deriving below our continuous model we 
assume that cells interact with one another subject to 
an excluded volume constraint. In the CPM a multidi- 
mensional integer index is associated with each lattice 
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site (pixel) to indicate that a pixel belongs to a cell of 
particular type or medium. Each cell is represented by 
a cluster of pixels with the same index. Pixels evolve 
according to the classical Metropolis algorithm based on 
Boltzmann statistics, and the effective energy 

E = E Adhesion + Ep erimeter + Epi e ld- (1) 

Namely, if a proposed change in a lattice configuration 
results in energy change AE, it is accepted with proba- 
bility 

where T represents an effective boundary fluctuation am- 
plitude of model cells in units of energy. Since the cells' 
environment is highly viscous, cells move to minimize 
their total energy consistent with imposed constraints 
and boundary conditions. If a change of a randomly cho- 
sen pixels' index causes cell-cell overlap it is abandoned. 
Otherwise, the acceptance probability is calculated using 
the corresponding energy change. If the change attempt 
is accepted, this results in changing location of the center 
of mass and dimensions of the cell. 

In this paper we assume that each cell has a rectangu- 
lar shape, that it moves or changes its shape by adding 
or removing a row or column of pixels (see Figure [T]) 
and that cells come into direct contact with each other. 
They also interact with each other over long distances 
through producing diffusing chemicals and reacting to 
local chemical gradients (process called chemotaxis) . Al- 
though we model adhesion between cells and the extra- 
cellular matrix (ECM), we neglect cell-cell adhesion and 
take into account cell-cell interaction from excluded vol- 
ume constraint meaning that cells cannot occupy the 
same volume. Under these assumptions terms in the 
Hamiltonian (TT]) have the following forms. EAdhesion 
phenomenologically describes the net adhesion or re- 
pulsion between the cell surface and ECM and it is a 
product of the binding energy per unit length, Jcm, 
and the length of an interface between the cell bound- 
ary and ECM: EAdhesion = 2J cm (L x + L y ). E P 

erimeter 

defines an energy penalty function for dimensions of 
a cell deviating from the target values Lp x and Lp ■ 

Eperimeter = A X (L X — Lp^) + A y (L y — Lp y ) where X x 

and A y are constants. Cells can move up or down gradi- 
ents of both diffusible chemical signals (i.e., chemotaxis) 
and insoluble ECM molecules (i.e., haptotaxis) described 
by Epi e id = iic(r)L x L y , r = (x,y) where c(r) is a local 
concentration of particular species of signaling molecules 
in the extracellular space and /x is an effective chemical 
potential. 

Let P(r, L, t) denote the probability density for a rect- 
angular cell with its center of mass at r to have dimen- 
sions L = (L x , L y ) at time t. We use vectors ei j2 to indi- 
cate changes in x and y dimensions: ei = Ar(l, 0), e 2 — 



Ar(0, 1). Let us normalize the total probability to the 
number of cells: J P(r, L, t)drdL = N. 

Now assume that cells cannot occupy the same space. 
This implies that position r' and size L' of any neigh- 
boring cell should satisfy the following excluded volume 

conditions: \x - x'\ > \y - y'\ > 

A discrete stochastic model of the cell dynamics un- 
der these conditions is described by the following master 
equation 

2 

P(r, L, t + e 2 At) = £ { [- - ^(r - |e, , L + ee,-; r, L, t) 

-$,>(r + |e i7 L + eej-,r, L, t) - T;(r + |e J; L - ee.,; r, L, t) 

-T r (r - Uj, L - eej ; r, L, t)] P(r, L, t) 

+$ J y(r,L;r+ -e^L- ee J ,t)P(r + ~ej,L - ee,-,i) 

+$j, r (r,L;r- ^e j: L - ee,-, t)P(r - -ej, L - ee^, t) 

+T ; (r,L;r- |e j7 L + ee,-, t)P(r - |e 3 , L + ee,-, t) 

+T r (r ) L;r+|e J -,L + eej,t)P(r+|e j ,L + ee i ,t)},(3) 

where T/(r, L; r', L', t) and T r (r, L; r', L', t) denote prob- 
abilities of transitions from a cell of length L' and 
center of mass at r' to a cell of dimensions L and 
center of mass at r without taking into account ex- 
cluded volume principle. (Terms with TJ and T r in 
Eq. ([3]) correspond to the case of decreasing cell size 
|L| < |L'| which justifies the neglect of excluded vol- 
ume.) 3>j,i(r, L; r', L', t) and $_y- )r (r, L;r',L',t) are prob- 
abilities of transitions taking into account excluded vol- 
ume. Subscripts / and r correspond to transitions by 
addition/removal of a row/colomn of pixels from the 
rear/lower and front /upper ends of a cell respectively. 
According to the CPM we have that T t (x, L; r', L') = 

T r (r,L;r',L') = |$(p(r, L) - E(r', L')) where the fac- 
tor of 1/8 is due to the fact that there are potentially 8 
possibilities for increasing or decreasing of L x and L y . 

We define ,(r, L; r', L') = T l[r) (r, L; r', L')[l - 
0j,r(/)( r ? where cj)j >r n\(r, L,i) is the probability of 

another cell being in the immediate neighborhood of a 
given cell and, therefore, preventing an increase of that 
cells' length or width (excluded volume). We neglect 
triple and higher order "collisions" between cells result- 
ing in the following approximation formulas 

0!, fc (r, L,t) - (iV-l)(eAr) 4 

L',a' 2 

2 , fc (r,L,<) = (A-l)(eAr) 4 

x J2Q(^±^-\x'-x\)P(r',L',t) Ly+L , (4) 

-w i I y — y+ s — 2 — 

L' ,x' 



where s — 1 for k = I, s = — 1 for k = r and factor 
is due to pairwise cell collisions. 

We found by using Monte Carlo simulations 
shown) that solutions of the master equation (E □ 
with general initial conditions quickly conver, 
P(r, L, t) = P Bo i tz (r,L)p(r,t) where P Bo itz(r,] 
Z{r)~ 1 exp(— pAEi eng th) is the Boltzmann dis 
tion and AE length = E(r,L) - E min = X x 
X y L 2 y + L x L v hc{y) and L = L L^ mm \ 
E min = i?(r, l/ mm )) is the minimal value c 
Hamiltonian JT]) achieved at L = LA™™) and Z 
(2eAr) 2 Eexp(-/3A^ en ^) ~ - ^ ■ 

^ 3 /3^/4A x A y -/i 2 c(r) 2 

is an asymptotic formula for a partition function. 

The typical fluctuation of cell dimensions (L x , L 
determined by /3X X ^L 2 ^ ~ 1. We now assume : 
dition that (5xq\ x ^> 1 and flUo^y ^S> 1, where a; 
yo are typical scales of P with respect to a; and y. mis 
means that xq^> L x , yo L y . We also assume that the 
concentration of chemoattractant c(r) is a slowly vary- 
ing function of r on a scale of the typical cell's length 
meaning that x c /L x 3> 1, Uc/L y 1, where x c and 
y c are typical scales for variation of c(r) in a; and y. We 
also make the additional biologically relevant assumption 
that 4X x X y ^> fi 2 c(r) 2 which means that change of typ- 
ical cell size due to chemotaxis (SL x chemo ' > , SLy Chemo ^ ) is 
small \5L^ c f h 1 mo h <§; Z^""™'. Under all above mentioned 
assumptions, the master Eq. (J3]) is transformed in the 
limit e -C 1 into an integro-differential equation describ- 
ing evolution of the probability density p(r, t) for the lo- 
cation of the cellular center of mass 

d t p = D 2 d 2 p - xod r ■ [p9rc(r)] 
+^-(N-l){d x [^ x p]+d y [i> y p]} 

= J W + L i mm) > f ') - p( x - L i mm) > 2/0] <V 

V, = | y + 4™")) - p(a/, y - L ( y m ^)]dx' 

I ) r> j (min) t (min) /r\ 

XQ = -D 2 fJ.pL x >L y ; ,(5) 

where D 2 = ^Xt^d 2 = d% + d%, X o = 
-D 2 ^L x mm) L y mm \ L x mm) = L Tsc - fy, 4™ m) = 
Lt v — ^j 2 ^ and J p(r)dr — N. Lastly, we couple this 
equation to an equation describing evolution of the ex- 
ternal (chemotactic) field c 

d t c = D c d 2 c — 7c + ap (6) 

where D c ,"f and a are diffusion, decay and production 
rates of the field respectively. Note that the chemical is 
produced by cells. 
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FIG. 2: Comparison between mesoscopic CPM and macro- 
scopic continuous model. (a) Plot of a two-dimensional 
probability density distributions for a CPM simulation of 12 
cells with e = 0.01 and numerical solution p(x, y, t) of the 
continuous Eq.©. (b) Cross sections of p cp m(xo,y,t) and 
Pcon{xo,y,t) at xo = 53.0 as functions of y. 



If excluded volume is not taken into account (i.e. as- 
suming ip x = i/jy = 0) Eqs. tp and ([6]) reduce to the 
classical Keller-Segel system [lj which has a finite time 
singularity and which was used for modeling collapse (ag- 
gregation) of bacterial colonies [2l[ . Addition of excluded 
volume significantly slows down collapse and, therefore, 
Eqs. iJSJ and ((6|) can be used for simulating cellular ag- 
gregation for a much longer period of time. Spongy bone 
formation considered in this paper, is accompanied by 
secretion of a viscous or solid ECM (see below) which 
quickly stabilizes a transient or mctastablc arrangement 
of cells into a persistent microanatomy and therefore also 
prevents collapse. 

Figure [2] demonstrates a very good agreement between 
typical CPM simulation and numerical solution of the 
continuous model ||5j) and ([6]). Both simulations were 
performed on a rectangular domain < x, y < 100 with 
simulation time t enc i — 100. Parameters were chosen as 
follows: Ar = 1, L Tx = L Ty = 3, X x = X y = 1.5, J cm = 
2, /? = 15, fi = 0.1, D c = 3.0, 7 = 0.00025 and a = 
0.2. The time interval between successive Monte Carlo 
steps was St = e 2 At = 0.0001, e = 0.01. Discrete form 
of the equation ^ was used to calculate the chemical 
field dynamics on a 200 x 200 lattice with the time step 
At c = 0.0125 and initial chemical field chosen in the 
form of c (x,y) = ( x ~ 70 ) ^(v- 60 ) _ The typical size of 
the mesh used in the continuous model was 1000 x 1000 
and the time step was 0.002. A large number of CPM 
simulations have been run to guarantee a representative 
statistical ensemble. We assumed that at each time step 
each cell released chemical content aAt c which was then 
distributed to four nearest chemical lattice sites. 

In what follows, we illustrate the efficacy of the model 
by applying it to the formation of spongy bone via 
the intramembranous route. In this developmental phe- 
nomenon, which generates portions of the skull, maxilla 
and mandible in vertebrate organisms, bone cells, or os- 
teoblasts, differentiate directly from loosely packed mes- 
enchymal cells. The differentiating cells secrete TGF- 
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FIG. 3: Simulation of spongy bone formation process. Ar 



1) L T:c = L Ty 
fi = -0.1, D c 



0.6, \ x = Ay = 1.5, J cm = 0.002, = 15, 
0.5, 7 = 0.014, At c = e 2 At = 0.01, e = 0.1, 
tend = 180. (a) Monte Carlo CPM simulation, a = 0.7. N = 
15000 cells were randomly distributed in a domain < x, y < 
100 with initial chemical field at zero, (b) Numerical solution 
of the continuous model resulting from a uniform initial cell 
density distribution and with 5% random fluctuation, a — 0.2. 
(c) Histological section of developing spongy bone in the rat 
skull. Trichrome stain. Photographed from a section in the 
New York Medical College Histology slide collection. The 
effective magnification of this image is about 2x that of a and 
b. Scale bar: 100 micrometers. 



beta which acts chemotactically, influencing cell migra- 
tion while simultaneously inducing production of ECM 
|22|, which in developing bone is termed osteoid 

Depending on local conditions, including initial cell 
density, the bone will progress to a dense state or stop 
at a spongy state, in which bony rods or trabeculae form 
a swiss-cheese-like network (see Figure 3c) that eventu- 
ally contains marrow tissue originating from the circu- 
lation. Our mesoscopic and macroscopic model simula- 
tions which start with initially dilute populations of cells 
in a chemotactic field, subject to an excluded volume 
constraint, result in a transiently appearing set of inter- 
connected multicellular trabeculae (see Figures 3a and 
3b) similar to the experimental picture (Figure 3c). In 
particular, in the simulations and the developing tissue 
there are many nodes from which three branches extend, 
but few with larger numbers. 

In summary, we have derived a macroscopic continu- 
ous model (O from a mesoscopic two-dimensional CPM 
with excluded volume constraint and coupled it to a 
model of chemotaxis ([To]). Numerical simulations confirm 
a very good agreement between the CPM and macro- 
scopic equations. Numerical analysis of the macroscopic 
model facilitated determination of conditions promoting 
formation of a lattice-like aggregation pattern. This per- 
mitted us to locate the parameter ranges within which 



the model cells in the CPM simulations behaved qualita- 
tively like the living cells that form multicellular branches 
in spongy bone by intramembranous ossification(Figure 
3c) . In contrast to earlier suggestions that the trabecular 
arrangement of spongy bone is based on pre-existing vas- 



cular patterns 23], or later- forming patterns of mineral 
deposition IjJ 20 j, our results suggest that it can arise 



from the self-organizing behavior of mesenchymal cells 
interacting with their ECM. 
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